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II 


ABSTRACT 


2.1 PREPARATION. An investigation was made as to pro- 
ducing a workable model for the transient analysis of a 
cryogenic hydrogen filling system. A series of programs and 
subprograms defining: the momentum, mass, and energy 
balances, the physical properties, the transport properties, 
and their interactions were devised. 


2.2 TEST. The program was modified for a simple theoretical 
test fluid. Exhaustive runs and modifications were made and 
at this point no stability has been achieved except in tri- 
vial cases. 
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Ill 


SUMMARY 


3.1 SCOPE AND PURPOSE. This investigation was of a 
theoretical nature. A pressure driven, flow controlled 
cryogenic filing system was modeled to accept time 
dependent input. The system consists of an upstream 
tank (feedtank), upstream piping, a control element, 
downstream piping, and a downstream tank (external tank). 
The piping configuration is contained in block data, while 
certain initial conditions, ambient temperature, run time f 
and time interval, are entered as input. 


3.2 RANGE OF VARIABLES. Though the algorithm is not limited 
in its variable range, some practical restrictions do exist. 
No inconsistency with the fundamental scientific constraints 
should be entered. Piping should be of order 1km, diameters 
0.2m to 0.5m, and node lengths lm to 10m. The time 
interval should be chosen such that it is no more than 10% 
of the residence time of the smallest node, however it 
should be long enough to avoid computational error. 

Pressures range from lE5Pa to 5E5Pa, while temperatures 
between 20K and 300K are considered. Mass flows of course 
should be consistent with the 10% of the residence time 
restriction. 


3.3 RESULTS. The transport and physical property rou- 
tines used previously were reprogrammed. Most of them 
worked fairly well while others did not work at all. 
Several were improved upon. A transient cryogenic 
hydrogen filling process was programmed and tested 
with a fictitious fluid. Except in trivial cases, the 
program failed to reach convergence as of this date. 
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Abbreviations and Acronyms List 


VARIABLE 

TYPE 

DEFINITION 

UNITS 

CP 

INTEGER*2 

NODE # OF CONTROL POINT 

LESS 

ET 

1! 

" " " EXTERNAL TANK 

II 

FRAME 

II 

TIME INTERVAL # 


K 

ft 

FITTING TYPE 

II 

NODE 

II 

NODE # 


PUMP 

ft 

NODE # OF PUMP 

II 

D 

REAL *8 

INSIDE DIAMETER 

M 

DELTAT 

•1 

TIME INTERVAL 

S 

DO 

If 

OUTSIDE DIAMETER 

M 

E 

tl 

U+v**2/2+G*Z 

M**2/S**2 

EPS 

ft 

PIPE ROUGHNESS 

M 

H 

II 

ENTHALPY 

M**2/S**2 

HO 

If 

H+v**2/2+G*Z 

M**2/S**2 

H 

ft 

HEAT TRANSFER COEFFICIENT 

KG/S**3/K 

F 

ft 

FANNING FRICTION FACTOR 

LESS 

K 

ft 

THERMAL CONDUCTIVITY 

KG*M/S**3/K 

KCO 

If 

TOTAL DISCHARGE COEFFICIENT 

LESS 

KELVIN 

tl 

BULK TEMPERATURE 

K 

KELWAL 

It 

WALL TEMPERATURE 

K 

M 

tl 

MASS OF NODE 

M 

MDOT ' 

tl 

MASS FLOW 

KG/S 

MU 

II 

VISCOSITY 

KG/M/S 

P 

fl 

PRESSURE 

KG/M/S* *2 

PF 

II 

FRICTIONAL PRESSURE DROP 

KG/M/S* *2 

PET 

It 

EXTERNAL TANK PRESSURE 

KG/M/S* *2 

PFT 

II 

FEED TANK PRESSURE 

KG/M/S* *2 

Q 

II 

HEAT FLUX 

KG/S * * 3 

QEXT 

It 

EXTERNAL HEAT FLUX 

KG/S* *3 

QMAX 

It 

MAXIMUM HEAT FLUX 

KG/S* *3 

REY 

II 

REYNOLDS NUMBER 

LESS 

RHO 

It 

DENSITY 

KG/M* *3 

U 

II 

INTERNAL ENERGY 

M**2/S**2 

TAMB 

II 

AMBIENT TEMPERATURE 

K 

TIME 

II 

ELAPSED TIME 

S 

TFINAL 

II 

TOTAL ELAPSED TIME 

S 

V 

II 

VELOCITY 

M/S 

VET 

II 

EXTERNAL TANK VOLUME 

M**3 

X 

II 

QUALITY 

LESS 

z 

II 

ELEVATION 
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VI 


BODY OF TEXT 


6.1 INTRODUCTION 

6.1.1 SUBJECT. The subject of this project is the transient 
study of a cryogenic line. The intricate procedure involved 
in the filling of the external tank on the shuttle, coupled 
the large capacity of a long filling system, suggests a 
complicated transient. Predicting the nature of this 
transient might enable NASA to develop a more efficient 
filling procedure. 

6.1.2 STATUS OF THE PROBLEM. NASA had developed some soft- 
ware to analyse transient flow of cryogenic fluids. However, 
all but one of these programs, fail to include the effects 
of heat transfer and two-phase flow. The program that in- 
cludes these effects the TCTP(l) or "transient cryogenic 
transfer program". TCTP has reasonable success analyzing 
LOX systems, but is not as reliable for LH2-H2 modeling. 

TCTP is not well structured and is very poorly documented, 
after careful analysis, it was decided to start anew rather 
than modify TCTP. 

6.1.3 SIGNIFICANCE OF WORK DONE. "H2FILL" , a highly struc- 
tured and rigidly documented FORTRAN code was developed to 
solve the problem. In the course of attempting to modify 
TCTP a set of subprograms defining the properties of 
hydrogen were developed. "H2FILL" is still being analyzed 
for run time stability at this point. 


6.2 MAIN TEXT-DESCRIPTIVE INFORMATION 

6.2.1 BACKGROUND. In the process of trying to decide 
whether to modify TCTP or start anew, the author in- 
investigated the possibility of locating software that 
might flowsheet the TCTP program and thereby facilitate 
the possibility of deciphering it. The "TAMU" code was 
located (2). It was available for $3500. After weighing 
the trade-offs between modification and ini tiation, i t 
was decided to go with the latter. 


6.2.2 CONFIGURATION. The configuration consists of the 
feedtank, the upstream piping, the control point (control 
valve), the downstream piping, and the external tank. 


The parameters of this configuration are described in the 
BLOCK DATA statement "H2BLOCK". The time sensitive para- 
meters are read in from an input list "H2LIS , LIS ' . 


6.2.3 PROGRAM. Program "H2FILL" the mainline program calls 
in four subroutines in sequence. These subroutines are: 

1 ) "H2INPUT" , 2)"H2INIT", 3 ) "H2M0MNRG" , AND 4 ) "H20UTPUT" . The 
data between the subroutines is transferred by an unlabeled 
COMMON. 


6. 2. 3.1 Input. Subprogram "H2 INPUT" reads in the initial 


inventories in the tanks, the ambient temperature, and 
the initial pressures in the tanks. Then, the time of the 
the run and the time interval are read. The remaining data 
read in are the tank pressures and the mass flow at each 
non-zero time interval. 


6. 2. 3. 2 Initiation. Subprogram "H2INIT" initializes the 


conditions in the nodes. "H2INIT" also establishes the init- 
ial properties through "H2LSAT" , "PROPHP",and "PROPPTGS" . In 
addition it returns elevation from the inventory through 
"FEEDTANK" and "EXTANK". 


6. 2. 3. 3 Processing. "H2M0MNRG" simultaneously solves the 


momentum, continuity, and energy equations; the properties 
and flows at each node are updated. "TRANS" updates the 
transport coefficients and the wall temperature at each node 
from the properties and the design data. 


6-. 2.3.4 Output. "H20UTPUT" creates a file "H2DAT.DAT" on 

which it prints the data in spreadsheets. The first set 
is a node scan of each time. The second set is a time scan 
of each node. 


6.3 MATHEMATICAL PRESENTATION 


6.3.1 MOMENTUM From Newton's second law (1) 
F=d ( M*V)/dTIME+ ( MDOT*V ) out- ( MDOT*V ) in 
where F is the net force. 
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The parameters of this configuration are described in the 
and 

M=RHO*PI/4*D**2*L (2) 

for a cylinder .Also note that 

dF«-PI/4*D**2* (dP+dPF+RHO*G*dZ ) ( 3 ) 

and 

MDOT«=RHO*V*PI/4*V*D**2 (4) 

combining (1), (2) (3), and (4) and integrating we get 
( P( NODE, FRAME- 1 )-P(NODE-l , FRAME-1 )+ 

( RHO ( NODE , FRAME-1 ) *V( NODE-1 , FRAME-1 ) * *2- 
RHO( NODE-1 , FRAME-1 ) *V( NODE-1 , FRAME-1 ) **2 )+ 

PF ( NODE , FRAME-1 ) + G* ( Z ( NODE ) -Z { NODE-1 ) ) =L ( NODE ) /DELTAT* 

( RHO ( NODE , FRAME-1 ) * ( V ( NODE , FRAME-1 ) - 
RHO (NODE, FRAME )*V( NODE, FRAME ) ) (5) 

Substituting (4) in (5) and rearranging the following 
appears : 

MDOT ( NODE , FRAME ) =MDOT ( NODE , FRAME-1 ) - 

PI/4 *D ( NODE ) * *2 *DELTAT/L ( NODE ) * 

( P ( NODE , FRAME-1 ) -P ( NODE-1 , FRAME-1 ) + 

RHO (NODE, FRAME-1 ) *V( NODE , FRAME-1 ) **2- 
RHO( NODE-1 , FRAME-1 ) *V( NODE-1 , FRAME-1 ) **2 
+PF ( NODE , FRAME-1 ) +RHO ( NODE , FRAME-1 ) * 

G* ( Z ( NODE ) -Z ( NODE-1 ) ) ) (6) 

which updates the mass flow. 


6.3.2 CONTINUITY. For mass conversation one sees that 
dM/dTIME=MDOTin~MDOTout (7) 

integrating we get 

M ( NODE , FRAME ) =M ( NODE , FRAME-1 ) +DELTAT* 

MDOT ( NODE- 1, FRAME) -MDOT (NODE, FRAME) ) (8) 

and 

RHO ( NODE , FRAME ) =M ( NODE , FRAME ) / 

(PI/4*D(NODE)**2*L(NODE) ) (9) 

also 

V ( NODE , FRAME ) =MDOT ( NODE , FRAME ) / 

( RHO ( NODE , FRAME ) *PI/4 *D ( NODE , FRAME ) * *2 ) (10) 

hence the momentum and mass conservation laws has 
allowed us to update MDOT, M, RHO, and V. However, in 
in order to ascertain the state of the system we need an- 
other independent variable. Energy conservation will satisfy 
this requirement. 


6.3.3 ENERGY. For the energy conservation (1) through (8) 
d ( MDOT* E)/dTIME=( MDOT* (U+V**2/2+G*Z ) ) in 

—MDOT* ( U+V* *2/2+G*Z) )out 

+ ( heat ) in- ( work ) out , ( 7 ) (11) 

where ( work ) out- ( MDOT*P/RHO ) out- ( MDOT*P/RHO ) in, ( 1 ) (12) 

and ( heat ) in=HC*PI *D*L* ( KELWAL- KELVIN ) , ( 1 ) (13) 
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combining (11), (12), and (13) and noting the definition 
of HO we obtain: 

d ( M*E ) /dTIME“ ( MDOT* HO ) in- ( MDOT*HO ) OUt+ 

HC*PI*D*L*( KELWAL— KELVIN ) , ( 7 ) . (14) 

integrating and solving for the update of E gives 

E ( NODE , FRAME ) = ( ( ( MDOT ( NODE-1 , FRAME-1 ) *H0 ( NODE-1 , FRAME-1 ) - 

MDOT ( NODE , FRAME ) *H0 ( NODE , FRAME ) + 

HC ( NODE , FRAME-1 ) *PI *D ( NODE ) *L ( NODE ) * 

(KELWAL (NODE, FRAME) -KELVIN (NODE, FRAME) ) )* 
DELTAT+M ( NODE , FRAME-1 ) *E ( NODE , FRAME-1 ) / 
M(NODE, FRAME) (15) 

then we update U 

U ( NODE , FRAME ) =E ( NODE , FRAME ) -V ( NODE , FRAME ) * * 2/2-G* Z ( NODE ) ( 1 6 ) 
With U and RHO determined the state of the system is also 
determined. Therefore P, H, KELVIN, MU, K and X can all be 
updated. A combination of physical properties an the flow 
conditions then can produce F, HC, and eventually KELWAL, 

PF, and HO. Hence the update would be complete, and the 
next iteration can proceed. 


6.4 DISCUSSION 


6.4.1 HISTORY. NASA has several programs that simulate 
the cryogenic fueling systems at the Kennedy Space 
Center. However, only one of those algorithms accounts 
for both heat transfer and two-phase flow. The program 
with both of these attributes is the TCTP code. TCTP 
seems to be satisfactory for oxygen, but fails in 
describing hydrogen fueling. The authors' suspicion is 
that the hydrogen properties are "stiff”, and are not 
converging with the numerical techniques used in TCTP. 


6.4.2 AGENDA. In order to produce a program that would 
successfully model the hydrogen fueling system, the question 
of whether to modify the existing TCTP, or to start anew had 
to be contemplated. Ordinarily, the easier path would seem 
to be to attempt modification of the existing code. However, 
TCTP turned out to be very poorly structured, and hardly 
documented at all. Due to the above, and the fact that some- 
thing might be gained from a new approach, it was decided 
that a new approach would be the way to go. 


6.4.3 PROJECT. A highly structured intricately documented 
fortran code was started and developed for numerically 
solving the momentum, continuity, and energy equations. 
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VII 


CONCLUDING SECTION 


7.X CONCLUDING REMARKS 


7.1.1 PHYSICAL PROPERTIES. A battery of subprograms 
describing the properties of cryogenic hydrogen were 
developed. The mode was made compatible to the SI system of 
units. Most of the data agreed reasonably well with NBS 
data, while did not agree well at all. On several points 
the author improved on the existing algorithms. Extensive 
documentation and perfect structure was used at all times. 


7.1.2 TRANSPORT PROPERTIES. An improved algorithm was prod- 
uced for flow conditions and the thermal equations were 
rewritten. 


7.1.3 TRANSIENT CRYOGENIC FLOW ALGORITHM. A highly 
structured, extensively documented program was written. A 
ficticious fluid was introduced and at this point conver_ 
gence to reasonable numbers has not as yet been achieved. 


7.1.4 RECOMMENDATIONS. The author that some minor revision 
would stabilize the run time condition of H2FILL. The 
author also believes that when actual cryogenic hydrogen 
properties are evaluated the previous troubles experienced 
with tctp will reappear. A more novel method of presenting 
the equation of state, might alleviate the problem. Differ- 
ential techniques with smoothing might provide the answer. 
Other alternatives in the numerical analysis of the equation 
of state might also be tried. One might also note that 
proper adjustment of the input parameters can also bring 
about convergence. 
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VIII 


APPENDIX 

8.8 APPENDIX A 


8.8.1 TIME SCAN OF NODES 


FEED TAJIK AT NODE « 0 

control point at node ■ 4 

EXTERNAL TANK AT MODE >35 
NOOE I - ) 


TIME 

0 . OOOOOOD+OO 

0.1500000-01 
0 . JOOOOOD-Ol 
0 . 4 50000D-0 1 
0 . < 000000-01 
0 . 7500000-0 1 
0 . 0000000-01 
0 . 1050000*00 
0.1700000*00 
0.1350000*00 
0.1500000*00 
0.1(50000*00 
0 . 1000000*00 
0.1050000*00 
0.3100000*00 
0. J25000D*00 
0 . 3400000*00 
0 . 355000 D* 0 0* 
0.2300000*00 
0.3050000*00 
n *nnfinnr\»(ti\ 


QUALITY 

0 . OOOOOOD+OO 
0 .OOOOOOD+OO 
0 . 0000000*00 
0 . 0000000*00 
0 . 0000000*00 
0 . 0000000*00 
0 . 0000000*00 
0 . 0000000*00 
0 . 0000000*00 
0.0000000*00 
0 . 0000000*00 
0 .000 0000*00 
0 . OOOOOOD+OO 
0 .OOOOOOD+OO 
0. OOOOOOD + OO 
0 .OOOOOOD+OO 
0 .OOOOOOD+OO 
0 .OOOOOOD+OO 
0 .OOOOOOD + OO 
0 .OOOOOOD+OO 
0 OOOOOOD+OO 


PRESSURE 

0.499421 D+00 
O'. 4994210*04 
0 . 499421D+O0 
0.4994 2 1D+06 
0. 499431 D* 00 
0.499421 D*06 
0.4994210+00 
0.499421D+0I 
0 .499421 D+0 4 
0.4994 2 1D+04 
0.499421D+04 
0.49941 70+04 
0.4994260*06 
0. 499696 D+06 
0.4996070*06 
0.412641 D+06 
0.453474D+06 
0 . 2641 7JD+07 
-0.126012D+OI 
0. 2395410+10 
-0.794325D+11 


TEMPERATURE 

0.203 967D+02 
0.2030670*02 
0.20306 7D+ 0 2 
0.2039670*02 
0.2039670+02 
0.2039670+02 
0 . 20 39670*02 
0.2039670+02 
0 . 2039670*02 
0. 203967 D+ 02 
0.2039670*02 
0.2039660+02 
0.2039690+02 
0.2040710+02 
0.2040420+02 
0 .1971450+02 
0 .1453760+02 
6 . 472791 D+ 02 
-0.467474D+03 
-0.66U94D+04 
0.432444D+07 


WALL TEMP. 

0.2039000+02 
0 .2039000+02 
0.203900 D+02 
0 .203900D+02 
0.2039000+02 
0 .2039000+02 
0.2039000+02 
0. 2039000+02 
0 . 2039000+02 
0. 2039000+02 
0.2039000+02 
0 . 20 3 9 000+02 
0.20 39000+02 
0.2039000+02 
0.2039000+02 
0.2039000+02 
0 .2039000+02 
0.2039000+02 
0. 2039000+02 
0 .2039000+02 
0 .2039000+02 


DENSITY 

0 , 54901SD+01 
0. 549015 D+ 01 
0 .5490150+01 
0 . 5490150+01 
0.3490150+01 
0.5190150+01 
0. 5490150+01 
0.5490150+01 
0.5490150+01 
0 . 5490150+01 
0.5490150+01 
0.5490150+01 
0.549015D+01 
0 .5490140+01 
0.54901 7 D+0 1 
0 . 5444030+01 
0 . 3444620+01 
0.6013070+01 
0.4444410+01 
-0 .4715010+00 
“0.4414100+01 


VELOCITY 

0. OOOOOOD+OO 
0.5404110-11 
0 .5404110-11 
0. 540411D-11 
0.5404110-11 
0. 5404110-11 
0.5404110-11 
0.5404110-11 
0.5404110-11 
0.5404110-11 
0.540 4 11 D-l 1 
0.5404110-11 
0.5404110-11 
0.540407D-11 
0.5404000-11 
0.5406050-11 
0.5409100-11 
0.5293630-11 
0.4901100-11 
-0.3652430-10 
-0.7211210-11 


HASS 

0.1150440+01 
0.1050440+01 
0.1050440+01 
0 . 1 6 50 4 40+01 
0.1050440+01 
0.1050440+01 
0. 1050440+01 
0 . It 504 40+01 
0.1150440*01 
0.1150440+01 
0 .1050440+01 
0 . 1I5044D+01 
0.1050440+01 
0.1950460+01 
0.1050450+01 
0 .144974D+01 
0. 1*41710+01 
0. 1*69060+01 
0.2037160+01 
-0. 2737900 + 00 
-0. 1306730+01 


MASS PLOW 

0. OOOOOOD+OO 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 100000 D -1 1 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0 . 1000000-11 
0.100000D-11 
0 . 100000 D -11 
0 . 1000000-11 
O.lOOOOOD-li 
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8.8.2 NODE SCAN OF TIME 


fit© TJLHJt At HOSE *1 6 

tOVTDOL >OJftT AT ftOJSI » 4 

tlTtlhil. AT HOPE - 33 

tUkMiL' TJ*C • ft.BB60eCD.O4 

■toot OOALITT llltiull ItRtimvIl WALL TIlCF. DUUTT trWCJTT KM I JUil »LOW 
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D . B D 0 D B 6 0.0 ft 
0 . fl 0 D D fi fl P* 0 D 

f . p&d ooeD.ce 
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0 . 0 D D D o D I>* ft ft 
ft. DOtlOODC.ftO 
0 . OCftO D 0l'« ft ft 
ft. DDftDOQb.ftO 
ft . 0 DOC 6 DXi.ftfl 
B . 0 0 6 D 0 OWftft 
6 . 0 D 0 DftO&.ftC 
0 . 0 0 6 0 o DO+ftft 
L . ODCOODt.flO 
fi.00QQ6CD.Cft 
ft .OP JvOC-JuOB 
fi.CC.fiQQDD.ftO 
0 .flDPPBPD.OB 
ft. OftDOQCCi.CiO 
ft .CDfifiOflb.Bfl 
B . DDOOOBC— ft® 
D . 0 D 0 0 0 OD.fl B 
B . 0 0 0 0 0 0 0*00 


• . )00 D fi 0 p.O 4 
d.umhd»P( 
fi . 4>942}D.G4 
ft. 49(42 ID. OB 
ft. J D 0 0 4 7 d*04 
0 . 0 Oft ft 4 1b«0 | 
B- JDBC47p.ft4 

• . JOftOt 7r>.ot 
0.2BB64TD.P4 
ft . JOftft ■ 70. CB 

• . J ftQft 4 Td» 0 t 

0 . ?flft»4 JD.D4 

• .2BBB4 TD.04 
B . 20P04TD-04 
ft.20ftft47D.0B 
B . 2966c)p.0| 
ft. )B ftft i 7 0.0 B 
B . 2ft 0 6 4 7 Q- 0 4 
C -20ft# 4 ?D.OB 
D . JOPPflp-.fif 
ft , 2flftO I 7D.04 
0 . iv D & 1 J P .0 i 
o.iftftftf 7&.OB 
0.209 9 47D.OB 
ft . JCftO I7D.0B 
9 . 20 P P 4 Tfi. 0 i 
ft . 2 B 0 ft® Td. 0 B 
ft . JCP04 7D.04 
ft . 2B8B4TD.0B 
D . 2 0 ft ft 4 7 D» 0 B 
0.20PB47D.04 
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8.9 APPENDIX B 


8.9.1 FORTRAN PROGRAM H2FILL 


program h2fill 

intege t * 2 cp,et,frame,kc<5,0:99) , node , pump 
real's cpr(0:99,0:99) , cvo ( 0 : 9 9 , 0 : 99 ) , 

4 <1(0:99) .delta t ,<10(0:99) ,#(0:99,0:99) ,eps(0:99) , 

4 h(0:99, 0:99), hc(0:99, 0:99), 

4 h0(0:99, 0:99), f(Q:99, 0:99), k(Q:99, 0:99), 

4 kco< 0:99, 0:99) , k el vin ( 0 : 9 9 , 0 : 99 ) , kel wal ( 0 : 99 , 

4 0:99), 1(0:99), a(0:99, 0:99) ,adot ( 0 : 99 , 0 : 99 ) , 

t, au (0:99, 0:99), p(0:99, 0:99), pf (0:99, 0:99), 

4 pet(0:99),pft(0:99),q(0:99,0:99), 

4 qext (0:99,0:99), 

4 qmax ( 0 : 99 ) , rey ( 0 : 99 ,0 : 99 ) , 

4 r ho (0:99, 0:99), taab , tf inal , time ( 0 : 99 ) , 

4 u(0:99, 0:99), v(0:99, 0:99), 

4 vet, x(0:99, 0:99), x(0:99) 

common 

4 cp,et ,frame, kc, node, pump, 

4 cpr.cvo, 

4 d , del tat , do , e , eps , 

4 h ,hc ,h0 , f ,k , 

4 kco , kelvin, kelual , 

4 l,a,adot, 

4 mu , p , pf , pet , 

4 pft ,q, qext , qmax, 

4 rey , rho , tamb , tf inal , 

4 time.u, 

4 v , vet , x , z 


H2FILL CALCULATES THE TRANSIENT OF A PRESSURE DRIVEN H2 FILL 
SYSTEM. THE PROGRAM REQUIRES THAT INITIAL CONDITION IN THE 
TWO TANKS, THE MASS FLOW AT THE CONTROL POINT, AND THE AMBIENT 
TEMPERATURE BE GIVEN. TIME INDEPENDENT DATA ARE INCLUDED IN A 
BLOCK DATA STATEMENT. A SINGLE UNLABLED COMMON IS USED IN ALL 
SUBPROGRAMS; FOR TRANSFER OF CONTROL. BOUNDARY CONDITIONS RE- 
QUIRED ARE' THE PRESSURES AT THE TOP OF THE TANKS , AND THE MASS 
FLOW AT EACH TIME INTERVAL. 


ALL DIMENSIONED VARIABLES ARE IN THE SI SYSTEM INTERNALLY. 

ON I/O NON-SI DATA ARE CONVERTED TO/FROM SI DATA AT THE I/O 
OR BLOCK DATA INTERFACE. 


VARIABLE TYPE 


DEFINITION 


CP 

ET 

FRAME 

KC 

NODE 

PUMP 

D 

DELTAT 

DO 

E 

EPS 

H 

HO 

HC 

F 

K 


INTEGER* 2 


REAI* 8 


NODE | OF CONTROL POINT 
" * * EXTERNAL TANK 

TIME INTERVAL « 

FITTING TYPE 
NODE | 

NODE « OF PUMP 
INSIDE DIAMETER 
TIME INTERVAL 
OUTSIDE DIAMETER 
U+V * * 2/2+G *2 
PIPE ROUGHNESS 
ENTHALPY 
H+V * * 2/2+G *Z 

HEAT TRANSFER COEFFICIENT 
FANNING FRICTION FACTOR 
THERMAL CONDUCTIVITY 
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c 

KCO 

n 

TOTAL DISCHARGE COEFFICIENT LESS 

* 

c 

KELVIN 

ft 

BULK TEMPERATURE 

K 

* 

c 

KELWAL 

« 

WALL TEMPERATURE 

K 

* 

c 

H 

ft 

MASS OF NODE 

M 

* 

c 

HDOT 

" 

MASS FLOW 

KG/S 

* 

c 

MU 

N 

VISCOSITY 

KG/M/S 

* 

c 

P 

M 

PRESSURE 

KG/M/S* *2 

* 

c 

PF 

" 

FRICTIONAL PRESSURE DROP 

KG/M/S* »2 

* 

c 

PET 

If 

EXTERNAL TANK PRESSURE 

KG/M/S * * 2 

* 

c 

PFT 

If 

FEED TANK PRESSURE 

KG/M/S* *2 

* 

c 

0 

n 

HEAT FLUX 

KG/S* *3 

* 

c 

QEXT 

H 

EXTERNAL HEAT FLUX 

KG/S* *3 

* 

c 

QMAX 

ft 

MAXIMUM HEAT FLUX 

KG/S* *3 

* 

c 

REY 

N 

REYNOLDS NUMBER 

LESS 

* 

c 

RHO 

" 

DENSITY 

KG/M* *3 

* 

c 

U 

" 

INTERNAL ENERGY 

M* * 2/S * * 2 

* 

c 

TAMB 

* 

AMBIENT TEMPERATURE 

K 

* 

c 

TIME 

n 

ELAPSED TIME 

S 

* 

c 

TFINAL 

« 

TOTAL ELAPSED TIME 

S 

* 

c 

V 

n 

VELOCITY 

M/S 

* 

c 

VET 

" 

EXTERNAL TANK VOLUME 

M**3 

* 

c 

X 

* 

QUALITY 

LESS 

* 

c 

z 

" 

ELEVATION 

M 

* 

c* 


A************************************************* 

c 

c* 

At**************************************************************** 

c 

INPUT REQUIRED: 



* 

c 



H2LIS.LIS 


* 

c 

OUTPUT GENERATED: 




c 



H 2 DAT . DAT 


* 

c 

SUBROUTINES 

CALLED: 



* 

c 



H2INPUT 


* 

c 



H2INIT 


* 

c 

* 


h 2MOMNRG 


* 

c 



H20UTPUT 


* 

c******************************************************************* 

c 

c* 

****************************************************************** 

c 

INFORMATION 

CONCERNING THE PIPING CONFIGURATION IS CONTAINED 

* 

c 

IN THE BLOCK DATA STATEMENT "H2BLOCK . BLK" . THE 

INITIAL AND 

* 

c 

BOUNDARY CONDITIONS 

ARE READ IN ON DEVICE #20 

FROM FILE 

* 


C "M2 LIST . LIS M . * 

CHH**MM»iiH«<i*****»HH**H*t)****M**********************il****** 

c 

£******************************************************************* 
C H2 INPUT IS THE INPUT SUBROUTINE. * 

£********************************«********************************** 
call h2input 

£******************************************************************* 
C H2INIT IS THE INITIALIZATION SUBROUTINE * 

^ft****A*(ifi*ft***ft****A*fi*A*ft******«*ft*A*********«**. ***************** 

call h2init 

£******************************************************************* 
C h2MOMNRG IS THE PROCESSING SUBROUTINE * 

C******************************************************************* 


call h2aoanrg 

C******************************************************************* 

C H20UTPUT IS THE OUTPUT SUBROUTINE * 

C************^****************************************************** 

call h2output 

£*****************************«************************************* 
C H20UTPUT CREATES A REPORT ON A NEW FILE "H2DAT.DAT" ON DE- * 
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nnnnnoonn nnnnnnnnnnnn n n n 


C VICE #21. 

C**********#***a***a**********************************#***********a 

C 

C *»******»«********»*M**«***M*********ft*************M«********i* 

C THEN 

c ». *.*.**. ...**..*..**.<...*..*************»»**«***»************»*»* 

stop 


AND 

•on*****************************************.*********.*********** 

and 

subroutine h2input 


K2INPUT READS "H2LIS.LIS". THE FIRST RECORD CONSISTS Or THE 
MASS IN THE FEED TANK, THE AMBIENT TEMPERATURE , THE INITIAL 
PRESSURE IN THE FEED TANK, AND THE INITIAL PRESSURE IN THE 
EXTERNAL TANK. THE SECOND RECORD CONTAINS THE TIME INTERVAL 
AND THE THE TOTAL TIME OF THE RUN. ALL OF THE 
REMAINING RECORDS CONSIST OF THE PRESSURE AT THE FEED TANK 
THE MASS FLOW AT THE CONTROL POINT, AND THE PRESSURE AT THE 
EXTERNAL TANK AT ALL TIMES. THE THE RECORDS CONSIST OF THREE 
FLOATING POINT NUMBERS EACH SEPERATED BV A COMMA. THE 
INITIAL MASS FLOW WILL BE SET TO 0.0. DURING INITIATION. 


integer *2 cp , et , f rame , kc ( 5 , 0 : 99 ) , node, pump 
real *8 cpr ( 0 : 99 , 0 : 99 ) , cvo ( 0 : 99 , 0 : 99 ) , 

& d(0:99),deltat, do (0:99), e(0:99, 0:99) ,eps(0:99) , 

& h(0 :99 ,0:99) ,hc(0:99 ,0 :99) , 

& hO (0:99, 0:99), {(0:99, 0:99), k(0:99, 0:99), 

& kco (0:99, 0:99) , kel vin ( 0 : 99 , 0 : 99 ) , kelwal ( 0 : 99 , 

8 0:99) ,1(0:99) ,m (0:99, 0:99) ,mdot (0:99, 0:99), 

8 mu (0:99, 0 : 99), p(0: 99, 0:99), pf (0:99, 0:99), 

& pet (0 :99 ) ,pft (0:99) , q ( 0 : 99 , 0 : 99 ) , 

t qext (0:99,0:99), 

8 qmax ( 0 : 99 ) , rey ( 0 :99 ,0:99 ) , 

8 r ho (0:99, 0:99), tamb , tf inal , time (0:99), 

8 u(0:99, 0:99) ,v(0:99, 0:99) , 

8 vet,x (0:99,0 :99 ),z (0:99) 

common 

8 cp, et, frame, kc, node, pump, 

8 cpr, cvo, 

8 d, deltat , do , e , eps , 

8 h , he , hO , f , k , 

8 kco , kel vin , kelwal , 

8 l,m,mdot, 

8 mu , p , pf , pet , 

8 pft,q, qext, qmax, 

8 rey , rho , tamb , tf inal , 

8 time,u, 

8 v,vet,x,z 


OPEN "H2LIS . LIS " AND CONNECT TO DEVICE I 20 * 

open ( f ile*'h21is . lis ' , unit»20, status*' old ' ) 

.ft**************************************************************** 

READ THE FIRST RECORD (INITIAL CONDITIONS) » 

read(20, *)»(!, 0), tamb, pft(0),pet(0) 


READ THE SECOND RECORD (TIME PARAMETERS) * 

A************************************************************ 

read(20,*)deltat,tfinal 
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£*** ftft **** ftft * ft ** ftft **** ft ** ftftft *** 6*1 

C READ THE REMAINING RECORDS (BOUNDARY CONDITIONS) * 

C**ft**A*ftftftft*ftfi**ft*«**ft*ft*4**A***ft*ft******************************* 


Cl 

& 

£ A A A A A 

c 

c * * * * * 
£ A A A A * 

C 

£ A A A A A 
£* A A A A 

c 

£ A A A A A 


read { 2 0 , * ) ( pft ( frame ) , Bdot ( cp , f ran* ) , 
pet ( f rame ) , 

f rame = l , tf inal/del tat ) 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA4 
DISCONNECT AND CLOSE THE FILE 

AAAAAAAAAAAAAAAAAAAAAAA&AAAAAi 


AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


lAAAAAAAA 


closa(unit«20,status='ke«p' 

*********************** 

THEN 

a********************** 
return 

************************* 

AND 

******************** 
end 

subroutine h2init 

£************************************** 

C - SUBROUTINE "H2INIT.FOR" SETS THE INITIAL CONDITIONS FOR * 

C H2 FILL . THE INITIAL CONDITION IS CONSIDERED TO BE STATIC. * 

C THE PRESSURE IS DETERMINED HYDROSTATICALLY, TAKING INTO * 

C CONSIDERATION THAT HO IS CONSERVED IN THIS INSTANCE, H IS * 

C DETERMINED, FORCING THE DETERMINATION OF THE OTHER PROP- * 

C ERTIES AT THE NODE. THE FEED TANK IS CONSIDERED TO BE * 

C PARTIALLY FILLED WITH SATURATED LIQUID, WHILE THE EXTERNAL * 

C TANK IS GAS AT AMBIENT TEMPERATURE. » 

C****************************************************************** 


c ********************************* 

C SUBROUTINES CALLED: 

C . FEEDTANK 

C H2LSAT 

C PROPHP 

C PROPPTGS 

£****************************************, 


( ********* 


integer*2 cp,et , frame ,kc(5,0:99) ,node , puap 
real*8 cpr(0:99, 0:99), cvo(0:99, 0:99), 
a d(0:99), do ltat, do (0:99), #(0:99, 0:99) ,#ps (0:99 ) , 

a h(0:99, 0:99), hc(0:99, 0:99), 

a h0(0:99, 0:99), f(0:99, 0:99), k(0:99, 0:99), 

a kco(0:99, 0:99), kelvin(0:99, 0:99), kelwal(0;99, 

& 0:99), 1(0:99), 8(0:99,0:99) , mdot ( 0 : 99 , 0 : 99 ) , 

a mu (0:99,0: 99), p(0:99, 0:99), pf(0:99, 0:99), 

a pet (0:99), pft(0:99),q(0:99, 0:99), 

a ge x t ( 0 : 9 9 ,0 :99 ) , 

a qmax(0:99),rey(0:99,0:99), 

a rho(0:99, 0:99), ta mb, tfinal, tine (0:99), 

a u( 0 : 99 , 0 : 99 ) , v ( 0 : 99 , 0 : 99 ) , 

a vet ,x CO : 99 ,0 :99 ) ,z ( 0 : 99 ) 

common 


a cp , et , t rame , kc , node , pump , 

a cpr,cvo, 

a d , deltat , do , e , ops , 

a h , he , hO , f , k , 

a kco , ke 1 vin , kelwal , 

a 1 , k , mdot , 

a mu , p ,p£ , pet , 

a pft , q , qex t , qaax , 

a ray , rho , tamb , tf inal , 

a time , u , 
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r 


t v , vet , x , z 

Q * *********************************************************** ****** 
C SIX VARIABLES ARE INITIALIZED * 

c * ************* **************************************************** 
p<0,0)«pft<0) 
p< et , 0 ) =pet ( 0 ) 
kelvin(0,0)=20.39d0 
kel vin ( et , 0 ) =tamb 
tine ( 0 )=0 . 
n ( 0 , 0 )=n ( 1 , 0 ) 
f rane=0 
nodesO 

c ************ ****************************************************** 
C "H2LSAT" RETURNS THE PROPERTIES OF THE SATURATED LIQUID * 

c * ***************************************************************** 
call h21sat 

c *, ,**,*. ***,.***************************************************** 
C "FEEDTANK" RETURNS THE ELEVATION OF THE LIQUID LEVEL * 

c * ********************************* ******************************** 
call feedtank 

c ****************************************************************** 
C FOR THE INITIAL CONDITION HO IS INVAR I ENT AND IT IS * 

C CALCULATED FOR THE TOP OF THE LIQUID LEVEL. * 

C * ***************************************************************** 
h0(0,0)=h(0,0)+9.81d0*z(0) 

c * ************** *************************************************** 
C FROM THE OUTLET OF THE FEED TANK TO THE INLET OF THE » 

C CONTROL POINT INITIAL PROPERTIES ARE CALCULATED. * 

c****************************************************************** 

do node=l,cp-l 

c * **************************************** ******************* ****** 
C THE PRESSURE IS CALCULATED HYDROSTATICALLY. * 

Q****************************************************************** 

p (node , p ) =rho ( node -1 , 0 ) 
i *9 . 81d0 *( z ( node ) -z (node-1 ) ) 

i +p(node-l,0) 

c * ******************************************************* ********** 
C H IS DETERMINED FROM HO. * 

c * ***************************************************************** 

h0(node,0)=h0(0,0) 

h(node,0)=h0(node,0)-9.81d0*z(node) 

£******•**********)*********************•********************»****** 

C " PROPHP " DETERMINES THE PROPERTIES FROM H AND P * 

£********************************************«**********»**»******* 
call prophp 

Q* *************** ******************************»*****»***********,* 

C FIVE MORE VARIABLES ARE INITIALIZED * 

c * ********************************** ******************************* 
ndot ( node ,0 ) = 0 . 

e(node,0)=h0(node,0)-p(node,0)/rho(node,0) 
he ( node , 0 ) =0 . 
pf ( node , 0 1=0 . 

kelwal(nod«,0)=kelvin(node,0) 
end do 

Q****************************************************************** 

C THE REMAINING MASSES OF THE UPSTREAM NODES ARE ESTABLISHED. * 

c , *************** ****** ********.*************, ******** ************* 

do node=2,cp-l 

■(node,0)=rho(node,0)*datan(l.dO) 

& *d ( node ) * *2 *1 ( node ) 

end do 

c * ***********»****,<,**********», **********************, ************ 
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C THE INITIAL CONDITION IS NOW ESTABLISHED PROM THE CONTROL 

C POINT TO THE EXTERNAL TANK. 

£**.***.*.**.****.****.. ***...* 

node=et 

^*Aid***«(i*il4****4**i(iiHHl*ftAiH 

C "PROPPTGS" RETURNS THE PROPERTIES OF A GAS GIVEN PRESSURE * 

C TEMPERATURE INPUT. » 

£*.****.*.**.. **A.*.*****.**ft***ft***.*.******ft********««***ft*.*.**. 

kelvin(node,0 )s>f aab 
call h21sat 

....ft.*.*..***.**.**.****.*****.***. ft****.**.***.*******.**** 

THE MASS AND HO FOR THE EXTERNAL TANK ARE ESTABLISHED. * 

...... ..... .ft....... .... **«**. ...ft***.*.*.*****.***********.. 

n( ot , 0 ) = rho (at , 0 ) *vet 

........ *.***. ***.***.. *..**. ...ft***.**.****.****** ********** 

FOR THE INITIAL CONDITION HO IS INVARIENT AND IT IS * 

CALCULATED FOR THE TOP OF THE EXTERNAL TANK. * 

************************************************************* 

h0(»t,0)=h(et,0)+9.81d0*r(et) 

************************************************************* 
FROM THE INLET OF THE EXTERNAL TANK TO THE OUTLET OF THE * 
CONTROL POINT THE INITIAL PROPERTIES ARE CALCULATED. * 

************************************************************* 

do node=et-l , cp ,-i 


£ * A * ft A 

C 

Q * * * * * 
Q+ * * * ft 

c 
c 

£ ft ft ft ft ft 
£ A ft ft ft ft 

c 
c 

£ ft <l A * ft 
£ * ft * ft A 

c 

£ ft * ft ft ft 


THE PRESSURE IS CALCULATED HYDROSTATICALLY. 

*ftftft****ftftftft*ft*ftA*ftftfi****ftftft**ftft***ft**< 


r ********** 


p ( node , 0 )«rho ( nod« + l , 0 ) 

& * 9 . 81d0 * ( z ( nod«+l )-z { node ) ) 

fc +p {node+1 , 0 ) 

(2ft**ft*ft*ft***********AAfc*ft*Aft* ************************************** 

C H IS DETERMINED FROM HO. * 

£*.**.. ft A*.*****.. * . A . A * . ********* * * * * ft.*..***** A***. .......ft ...... 

h0(n4de,0)=h0(et,0) 

h(node,0)=h0(noda,0)-9.81d0*s(noda) 

0 ****. ft*...*.**..*.*..*.....****.**********. *.*.*. *************. ... 

C "PROPHP" DETERMINES THE PROPERTIES FROM H AND P * 

call prophp 

ft. ft. ..ft*..*.*.*..*.*****.***********. ...ft**.. .****. .* 

C FIVE MORE VARIABLES ARE INITIALIZED * 

C*. ....A...... 

mdot ( node , 0 ) =0 . 

e (node , 0 )=h0 (node , 0 )-p (nod# , 0 )/rho (node , 0 ) 
he ( node , 0 ) =0 . 
p£ ( node , 0 ) .0 . 

kelwal ( node ,0)=kelvin(node,0) 
end do 

C* *.****... 

C THE REMAINING MASSES OF THE DOWNSTREAM NODES ARE 

C ESTABLISHED. 

£.**..***.. ft*. ....ft... ..********.. 


do node=cp+l , et-1 

ffi(nodo,0)-=rho(nod«,0)' , 4.d0*d*tan(l.d0) 
t *d ( node ) *1 ( node ) 

end do 

do frame»l , tf inal/del tat 
p ( 0 , frame )=p£ t ( f raae ) 
nodesO 

call h21sat 

p(*t, frame )"=pat( £ raae ) 

nod«««t 
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call h21sat 
and do 

do nod««0 ,at 

■dot ( noda , 0 )=mdot ( node , 1 ) 
end do 


THEN 


return 


AND 


and 

subroutine h2aomnrg 


" H 2M0MNRG " BASICALLY SOLVES THE TRANSIENT MOMENTUM AND 
ENERGY EQUATIONS IN ORDER TO UPDATE M, MDOT , RHO , V, E, U, 
AND THEN THE PHYSICAL PROPERTIES. AFTER THE ABOVE HAS 
BEEN ESTABLISHED, THE TRANSPORT COEFFICIENTS HC, F, AND KCO 
ARE DETERMINED. 


SUBROUTINES CALLED: 

EXTANK 
FEEDTANK 
PROPURHO 
TRANS 

***************************************************************** 
integer*2 cp, at, frame, kc(5, 0:99) , noda , pump 
real * 8 cpr(0:99,0:99) , cvo ( 0 : 9 9 , 0 : 99 ) , 
k d ( 0 : 9 9 ) ,deltat,do(0:99) ,e(0:99,0:99) ,eps(0:99) , 

k h(0:99, 0:99) ,hc(0:99, 0:99) , 

k hO (’0:99, 0:99), f(0:99, 0:99), k(0: 99, 0:99), 

k kco (0 :99,0:99) ,kelvin(0:99,0 :99 ),kelwal(0:99, 

k 0:99), 1(0:99), >(0:99,0:99) , Bdot ( 0 : 99 , 0 : 99 ) , 

k au (0:99, 0:99), p(0:99, 0:99), pf (0:99, 0:99), 

k pet ( 0 : 99 ) ,pft (0:99 ) , q ( 0 : 9 9 , 0 : 99 ) , 

k qext ( 0 : 99 , 0 : 99 ) , 

k qmax ( 0 : 99 ) , rey ( 0 : 99 , 0 : 99 ) , 

k rho (0:99,0:99) , tanb ,tfinal, tine (0:99) , 

k u(0 :99,0 :99 ) , v(0 :99 ,0 : 99 ) , 

k vet, x<0 :99, 0 :99) ,z(0:99) 

common 

k cp , at , f rase , kc , node , pump , 

k cpr.cvo, 

k d, deltat , do ,e , eps , 

k h ,hc , hO , f , k , 

k kco , kel vin , kelwal , 

k l,ra,mdot, 

k mu , p , pf , pet , 

k pf t , q , qext , qmax , 

k rey , rho , tamb, tf inal , 

k time,u, 

k v , vet , x , z 


START THE CLOCK. * 

do f rame=l , tf inal/deltat 


RECORD THE TIME. * 
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time(frane)*frais 0 *deltat 

C***ft***ft*ft*ftftftftft**ft***ft!ll*ft*6*****«******************************** 

C HO ACCUMULATION OCCURS WITHIN THE CONTROL POINT NODE. * 

C****************************************************************** 

adot ( cp-1 , frame )=adot ( cp , f raae ) 

C THE FEED TANK PRESSURE IS ZEROED IN ORDER TO ENTER THE DO * 

C LOOP. IT IS LATER UPDATED TO WITHIN 1% Or PFT IN AN * 

C ITERATIVE PROCESS. * 

C********************************************************'********** 

p ( cp-1 , f raae )«p ( 0 , f raae ) * 9 . d-1 
C do while 

C t. (( aba (p(cp-I , f raae )/p(cp-l , fraae-1 ) ) 

C i .It. . 99d0 ) .or. 

C & ( abs ( p < cp-I , f raae )/p< cp-1 , frame-1 ) ) 

C t .gt. l.OldO)) 

£******ft****ft********A********************************************* 

C FOR THE INTERMEDIATE PIPING UPSTREAM OP THE CONTROL POINT. * 

C************************************************************* * * * * * 

do node-cp-2 ,1,-1 

^*#*A4A*»ft4Jift(i*<!4***4(lii*49ft<i*4*<li****A****ft****<l*<l*lMS*<i*»**lHl***** 

C CALCULATE THE CROSS SECTIONAL AREA. * 

£******************66********************************************** 

s=datan(l.d0)’ > d(nod9 + l)**2 

£*********************************** ********************* ********** 
C UPDATE MDOT * 

£****************************************************************** 
m dot ( node, fram«)=adot (node, f raae-1 )-s *del tat/1 ( node+1 ) * 
t (p(node+I, frame-1) -p (node, frame-l)+ 

& v ( node+1 , f raae-1 )* *2* 

i rho ( node+1 , f raae-1 )- 

fc , v ( node , f raae-1 )** 2 * rho ( node , 

k t raae-1 ) +rho ( node+1 , f raae-1 ) 

i *9 . 8 ldO* ( z < node + 1 )-z (node ) )+pf ( 

k node+1 , fraoe-1 ) ) 

^************************* ***************************************** 
C UPDATE M, RHO, AND V » 

C*****i*****ft*rt*ftt»*<!**********ib*******ft***********«k****************ft 

end do 

do node«=cp-l , 2 , -1 
s=datan(l.d0)”d(node)**2 

® ( node , f ram a ) = ( ado t ( node-1 , f raae-1 ) — sdot ( node , f rane-1 ) ) 

I *de 1 tat+a ( node , f raae-1 ) 

f ho ( node , f cane )=m (nod« , frame )/s/l (node ) 
v ( node , f raao ) =mdot ( node , f rame )/s/rho ( node , f raae ) 
end do 

C*******ftftft*ft»*»ft**ft**A*ft************************ri 

C UPDATE MDOT FOR THE FEED TANK 

C******ft*ft**ft*fi****ftft*A**ae**ft*ft*******Aft********ri 

mdot ( Q« fra raa ) =mdot ( 1 , frame) 

a(0, frame) =m(0, f raae-1 )-adot ( 0 , f raae-1 ) * del tat 
£****************************************************** 


C UPDATE M, RHO, AND V FOR THE TERMINAL UPSTREAM NODES * 

C******************************************************************* 

m(l, frame) «m(0, frame) 

C.... THIS IS A COURSE APPROXIMATION FOR AN IDEAL GAS ONLY 

C A BARAMETRIC SUBROUTINE WOULD BE DEVELOPED FOR THE 

C . . . . GENERAL CASE 

rho(l, frame )>*rho(2, frame) 
v(l,frama)»v(2, frame) 

C******************************************************************* 

C UPDATE E,U,H,H0, AND PF FOR ALL THE UPSTREAM PIPE NODES * 
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do node=cp-l , 1 ,-l 

if (nodi .no. l)then 

« (nodi , f [«>« )>( ( adot (node-1 ,fraae-l)*hO( node-1 , 
f raae-1 )-adot (no da , f raae-1 ) *h0 
(node , f raae-1 ) )+hc (node , f raae-1 ) 

* ( 4 . dO * da tan ( 1 . dO ) *d (node ) »1 ( node ) ) » 
( kelwal (node , f raae-1 )-k a Ivin (node , 
f raae-1 ) ) *deltat+a( node , f raae-1 ) * 
e(node,fraae-l) ) /a (node, fraae) 

else 

e (node , fraae )»-adot (node , f raae-1 ) *h0 (node , f raae-1 ) 
+e (node , f raae-1 ) *a(node , f raae-1 )/ 
a(node , fraae ) 

endif 

u (node , fraae )»e (node , fraae ) 

-v(node,fraae)**2/2-9.81d0*z(node) 


"PROPURHO" RETURNS THE PHYSICAL PROPERTIES FROM AN 
INPUT OF U AND RHO . 


* 

* 


call propurho 

hO (node , fraae )»h (node , fraae )+v(node , fraae ) **2/2 
a +9 . 81d0*z (node ) 


"TRANS" RETURNS THE TRANSPORT PROPERTIES F , KCO , AND HC 
FROM THE PHYSICAL PROPERTIES AND CERTAIN FLOW 
PARAMETERS THAT HAVE ALREADY BEEN ESTABLISHED. TW IS 
ALSO RETURNED. 

i f ( node .ne. 1 ) then 
call trans 

pf (node , f rams )=■ ( 4 .dO *f (node , fraae )+kco (node , fraae ) ) * 
a rho ( node , f rame ) * 

a v (node , frame ) **2/2* 

a 1 ( node ) /d ( node ) 

endif 
end do 


"FEEDTANK" RETURNS THE ELEVATION OF THE LIQUID LEVEL IN THE * 
FEED TANK FROM THE FEED TANK MASS. * 

******************************************************************* 

call feedtank 


UPDATE E,H,AND HO FOR THE TOP OF THE FEED TANK 


e ( 0 , fraae )« (-adot ( 1 , f raae-1 ) "hO ( 1 , f raae-1 ) * 
a deltat+a( 1 , f raae-1 ) *e ( 0 , f raae-1 ) ) 

a /a(l, frame) 

u(0, fraae )>e(0, fra me )- 
a . 9.81d0*z(0) 

node«0 

call propurho 
hOtO.fraaelBhfO.fraae)* 
a 9 . 8 ldO * z ( 0 ) 


THIS IS 

THE THROTTLING 

CONDITION. 

* 

hO ( cp , 

, fraae ) »h0 ( cp-1 , 

, fraae ) 



INITALI ZE P INORDER TO ENTER LOOP. AN ITERATIVE PROCESS WILL* 
THEN BRING UPDATE IT TO WITHIN 1% OF PET. * 
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(;»***»***»*•***«»«*»*»*«»»»****»»*»»»*«*******************»*»****** 


p ( at , frame ) =0 . 
do While 


& ( aba ( pa t ( frame ) -p ( at , frame ) ) 

i .gt. 0 . 01d0*pet ( frame ) ) 

£************************************<**♦****************#********* 
C EVALUATE M DOT, RHO, V, M, E, HO, U, H, PLUS THE PHYSICAL 

C AND TRANSPORT PROPERTIES FOR THE INTERMEDIATE NODES. NOTE 

C THAT TW IS ALSO RETURNED. 

£************<*******«£>******«************************************* 


k 

k 

k 

k 

k 

k 

k 

k 


k 

k 

k 

k 

k 

k 

k 


k 

k 


do node»®t-l ,cp+l ,-l 
s=da t an (l.d0)*d(node)**2 

mdot(node,frame)>»adot(node, frame-1 )-s * da It at/1 (node ) * 
(p(nod®,frame-l)-p(nod®-l,fraae-l)+ 
v(node,£rame-l)**2* 
rho(node,fraae-l)- 
v(node-l,frarae-l)**2*rho( node-1 , 
frame-1 )+rho(node-l,fram®-l) 
*9.81d0*(z(node)-z(node-l))+pf( 
node , f rame-1 ) ) 

a ( node , frame ) ■»< ndot ( node-1 , f raae-1 )~adot (node , f rame-1 ) ) * del tat 
+m< node, frame- 1) 

rho( no de,frame)=m(node, frame ) /s/1 (node ) 
v ( node , frame )=radot (node , frame 9/s/rho (node , frame ) 
e(node,frame)=((adot(node-l,£rame-l)*hO (node-1 , 
frame-l)-mdot(node,frame-l)*hO 
( node , f rame-1 ) )+hc (node , frame-1 ) 

* ( 4 . dO * da tend. dO ) *d ( node ) *1 ( node ) ) * 

( kelwal (node , frame-1 ) -ke Ivin ( node , 
frame-1)) * del tat+a ( node , f rame-l ) * 
o (node, frame— 1) )/m( node , frame ) 
u( no de,frame)=e(node, frame) 

-v(node, frame )** 2/2-9. 81d0*z(node) 

call propurho 

h ( node , f rame )=*u ( node , frame )+p(node , f raae ) 

/rho ( node , f rame ) 

ho ( node , frame )=h (node , frame ) -tv (node, frame ) **2/2 

+9 .81d0*z(node) 


call trans 

pf( node, frame )=(4.d0*f(node, frame )+k co (node , frame ) ) * 
t rho ( node , f rame ) * 

t v ( node , f rame )* *2/2 * 

a 1 ( node ) /d ( node ) 

end do 

C****ftft**ftft**ftft**ftft*ftftftft*ft*ft!h***ft********* **************** ********* 

C UPDATE THE DOWNSTREAM PROPERTIES FOR THE CONTROL ELEMENT * 

C AND THE TOP OF THE EXTERNAL TANK. * 

£****ft**il****ft*ft*Aaft*ft*ft*****ft************** ********* ************** 


rho ( cp , frame )= rho ( cp+1 , f rame ) 

v(cp, frame )=mdot(cp,frame)/rho(cp, frame) 

& /datan ( 1 . dO )/d ( node ) * * 2 

h(cp, fra me )=h0(cp, frame) 
t -v ( cp , f rams )* *2/2 

a -9 . 8 ldO *z ( f rame ) 

node»cp 

************************************************ ItltltlKllltilttll 

"PROPHRHO" RETURNS THE PHYSICAL PROPERTIES FROM AN * 

INPUT OF H AND RHO. * 

a*************************************** ********** *•«****••••****•• 

call prophrho 
e(cp, frame )»hO(cp,frame)~ 
t p ( cp , f rame )/ 
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k rho ( cp , f rame ) 

u(cp,frame)=h(cp, frame )- 
& p(cp, frame)/ 

& rho ( cp , frame ) 


"EXTANK" RETURNS THE ELEVATION OF THE FLUID LEVEL FROM THE 
MASS AND QUALITY IN THE EXTERNAL TANK. 


call extank 

e(et,frame)=ndot(et-l,frame-l)*hO(et,frame-l) 

6 * de 1 t a t +a ( e t , frame-1 ) *e ( et , f rame-1 ) 

k /m(et, frame) 

u(et, frame )=e(et, frame )- 
k 9 . 81d0 * x ( et ) 

nodeaet 
call propurho 
hO ( et , f rame ) =e ( et , f rame ) + 
k p(et, frame)/ 

k rho ( et , f rame ) 

h ( et , f rame ) *u ( et , frame ) + 

k p(et, frame)/ 

k rho ( et , f rame ) 

end do 

eft****************************************************************** 

C RESET THE CLOCK * 

£******•***********»*********************************< 
end do 

(;»•»»»*•»*«»*»**«*»**********»( 

C THEN 

C*»»**»**»**«**«*»****»»*****»*l 

return 

^*****•***********•************************•********** 4 ************* 

C AND * 

£***************&******************< 
end 

subroutine h2output 

C "H20UTPUT" CREATES THE FILE "H2DAT.DAT" ON DEVICE # 21. IT • 

C THEN PRINTS A TABULAR OUTPUT OF THE TIME AT VARIOUS NODES, * 

C FOLLOWED BY A REARRANGEMENT OF THE DATA FOR THE NODES AT * 

C VARIOUS TIMES. * 




integer*2 cp , et , frame , kc ( 5 , 0 : 99 ) , node , pump 
real * 8 cpr(0:99,0:99) , cvo ( 0 : 9 9 , 0 : 99 ) , 
k d(0 :99 i .delta t , do ( 0 :99 ) , # (0 :99 ,0:99 ) ,eps (0:99) , 

k h(0 :99 ,0 :99 ) ,hc(0 :99 ,0 :99 ) , 

k ho (0:99, 0:99), f(0 : 99, 0:99), k<0:99, 0:99), 

k kco ( 0 : 99 , 0 : 99 ) , k el vin ( 0 : 9 9 , 0 : 99 ) ,kelva 1(0:99, 

k 0:99), 1(0:99), *(0:99, 0:99) , mdot ( 0 : 9 9 , 0 : 99 ) , 

k mu (0:99, 0:99), p(0:99, 0:99), pf (0:99 , 0:99), 

k pet (0:99), pft(0:99),q<0:99, 0:99), 

k qext ('0 : 99 , 0 : 99 ) , 

k qmax (0:99), rey (0:99, 0:99), 

k rho(0:99,0:99),tamb,tfinal,time(0:99), 

k u(0:99, 0:99), v(0:99, 0:99), 

k vet, x(0 : 99, 0 :99 ) ,z(0:99) 

common 

k cp , et , f rame , kc , node , pump , 

t cpr,cvo, 

k d , deltat , do , e , eps , 

k h ,hc ,h0 , f , k , 

k kco , ke Ivin , ke lwal , 
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4 l,m,mdot, 

4 mu , p , pf , pa t , 

4 pft,q,qext, qnax , 

k rey , rho , tamb , t f inal , 

k tima.u, 

k v,v«t,x,i 

£*************•************************ 

C CREATE OUTPUT FILE "H2DAT.DAT" AND CONNECT TO DEVICE « 21. 


open ( f ile= ' h2dat . dat ' , unit=21, status® ' new ’ ) 
£****************************************************************** 
C START THE CLOCK. 


do frame»0 , tf inal/del tat 

C*************************************** 

C PRINT THE HEADING. 


*************** 


write ( 21 , 100 ) cp , at , tine ( f raaa ) 

100 format ( lx , 'FEED TANK AT NODE = 0',/, 

k lx, 'CONTROL POINT AT NODE = ’,i2,/, 

k lx, 'EXTERNAL TANK AT NODE « ',i2,/, 

4 lx, 'ELAPSED TIME ■ ',dl3.6,/> 

write { 21 , 101 ) 

101 f ormat ( lx , 5x , 'NODE ' , 5x , 

k 4x , ' QUALITY ' 3x , 

k 3x, 'PRESSURE' ,3x, 

4 2x, 'TEMPERATURE' , lx, 

4 2x , 'WALL TEMP . ' , 2x , 

4 4x, 'DENSITY' ,3x, 

4 3x, 'VELOCITY' ,3x, 

4 5x , 'MASS ' , 5x , 

4 3x , 'MASS FLOW' ,/) 

£***********************************< 

C SCAN THE NODES. 

C************************************* 

do node=Q,et 

write ( 21 , 102 ) noda , x ( node , frame ) ,p ( nod# .frame ) , 

4 kelvin(noda, frame), kelwal(node, frame), 

4 rho ( noda , f raae ), v ( node , f rane ) , 

4 m ( node , f raae |, mdot ( node , frame ) 

102 format(lx,6x,i2,6x,8(lx,dl3.6) ) 

C********************************************************** 

C NEXT NODE 


***************** 


and do 

C********************* 

C RESET THE CLOCK 


**************************************** 


and do 

************************************************! 
TOP OF THE PAGE FOR THE TIME SCAN OF THE NODES. 


write (21,103) 

103 format ( ' 1 ' ) 

C****‘ ************ **.*•****»< 

C SCAN THE NODES. 

£**********************************< 
do node=0,et 


PRINT A HEADING. 

************************ 6*4 

wr i t® ( 21 , 104 ) cp, at, node 
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104 f ormat ( lx , ' FEED TANK AT NODE a 0',/, 

k lx, 'CONTROL POINT AT NODE = 

k lx, 'EXTERNAL TANK AT NODE = ',i2,/, 

i lx, 'NODE » ■ ' ,i2,/) 

wri to ( 21 , 105 ) 

105 f ormat ( lx , Sx , 'TIME ', 5x , 

k 4x, 'QUALITY '3x, 

k 3x, 'PRESSURE' ,3x, 

4 2X , 'TEMPERATURE' , lx, 

k 2x , 'WALL TEMP. ' ,2x, 

k 4x , ' DENSITY' , 3x , 

k 3x, 'VELOCITY' ,3x, 

t 5x, 'MASS' ,5x, 

k 3x , 'MASS FLOW' ,/) 


C* 

C 


START THE CLOCK 


do £ rame=0 , tf inal/del tat 


PRINT THE DATA. 

writa(21,106)time( frame ) , x (node , f rame ) , p ( node , frame ) , 


k kelvintnode, frame), kelwal( node, frame), 

k rho ( node , f rame ), v ( node , f rame ) , 

i a ( node , f rame ) ,mdot ( node , f rame ) 

106 f orma t ( lx , 9 ( lx , dl 3 . 6 ) ) 


C»* 

C 

C * * 

c** 

c 

c*» 

c* * 
c 

C»* 


c 

c** 

c** 

c 

c** 


*************************************< 

RESET THE CLOCK. 
******************************1 

end do 

*************************4 

NEXT NODE. 

*****************************< 
end do 

**tk***ft********ft***ft*< 

CLOSE AND DISCONNECT FILE "H2DAT.DAT". 

**********ftftft**ft«*********<M 


close ( unit»21 ) 


return 


AND 


* * * » * 
* 

* * * * * 


k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 


end 

subroutine h21sat 

integer*2 cp, et, frame, kc(5, 0:99) , node, pump 
real *8 cpr(0:99,0:99) , cvo ( 0 : 9 9 , 0 : 9 9 ) , 

d(0:99), del tat, do (0:99), #(0:99, 0:99) ,eps (0:99) 

h(0: 99, 0:99), he (0:99, 0:99), 

hO (0:99, 0:99), f(0:99, 0:99), k(0:99, 0:99), 

kco <0:99, 0:99) , kel vin ( 0 : 99 , 0 : 9 9 ) , kelwal ( 0 : 99 , 

0:99), 1(0:99), m(0:99, 0:99) , mdot ( 0 : 99 , 0 : 99 ) , 

mu< 0 : 99, 0 : 99 ),p< 0 : 99,0 :99),pf (0:99, 0:99), 

pet(0:99),pft(0:99), <1(0:99, 0:99), 

qext(0:99,0:99 ) , 

qm&x (0:99) ,r«y (0:99, 0:99), 

rho (0:99, 0:99) ,tamb , tf inal , time ( 0 :99 ) , 

u(0 : 99,0: 99), v(0:99, 0:99), 

vet, x(0:99, 0:99) , z ( 0 : 99 ) 


common 
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a cp , et , f rame , kc .node , pump , 

& cpr.cvo, 

& d , deitat , do , • , aps , 

a h , he , hO , f , k , 

a kco , kelvin , kelwal , 

a 1 , a , mdo t , 

a au,p,pf,pet, 

a pf t ,q , qext ,q*ax , 

a ray , rho , tamb , tf inal , 

a time , u , 

a v , vet , x , z 

call propptgs 
return 
end 

subroutine prophrho 

integer*;! cp.et, frame, kc(5, 0:99) .node .pump 
r ea 1 * 8 cpr ( 0 : 9 9 , 0 : 9 9 ) ,cvo(0:99.0:99>, 
a d ( 0 : 9 9 ) , deitat , do ( 0: 99 ),e(0:99, 0:99), eps (0:99), 

a h ( 0 : 99 , 0 : 99 ) ,hc( 0 : 99 , 0 : 99 ) , 

a h0(0 -.99,0 : 99), f(0:99, 0:99), k(Q:99, 0:99), 

a kco ( 0 : 99 , 0 :99 ) , ke 1 vin ( 0 : 9 9 , 0 : 99 ) .kelwal ( 0 : 99 , 

a 0:99), 1(0:99), a(0:99, 0:99) , adot ( 0 : 99 , 0 : 99 ) , 

a mu (0:99, 0:99), p(0:99, 0:99), pf (0:99, 0:99), 

a pet (0:99) ,pft (0:99) ,q(Q:99, 0:99), 

a qext ( 0 : 99 , 0 : 99 ) , 

a qmax ( 0 : 99 ) , rey ( 0 : 99 , 0 : 99 ) , 

a rho(0:99,0:99),tamb,tfinal,time(0:99), 

a u ( 0 : 99 , 0 : 99 ) , v( 0 : 99 ,0 :99 ) , 

a vet, x(0:99, 0:99), 2(0:99) 

common 


a 

a 

a 

a 

a 

a 

a 

a 

a 

a 

a 


a 

a 

a 

a 

a 

a 

a 

a 

a 

a 

a 

a 


cp, et, frame, kc, node, pump, 
cpr , evo , 

d;deltat,do,e,eps, 
h , he ,h0 , f , k , 
kco, kelvin, kelwal, 

1 , m , mdot , 
mu ,p,pf , pet , 
pft,q, qext, qmax, 
rey , rho , tamb , tf inal , 
time , u , 
v , vet , x , z 

u(node, frame) =h (node, frame)/!. 4d0 

call propurho 

return 

end 

subroutine propptgs 

integer * 2 cp,et, frame, kc(S, 0:99), node, pump 
real *6 cp r ( 0 : 9 9 , 0 : 9 9 ) , c vo ( 0 : 9 9 , 0 : 9 9 ) , 

d(0: 99), deitat, do (0:99), e(0:99, 0:99) ,eps (0:99) 

h(0 -.99,0:99), he (0:99, 0:99), 

h0( 0-: 99, 0:99), f(0 : 99, 0:99), k(0:99, 0:99), 

kco(0 :99,0:99) , k e 1 vin ( 0 : 9 9 , 0 : 9 9 ) , kelwal ( 0 : 99 , 

0: 99), 1(0:99), m(0:99, 0:99) ,mdot (0:99,0:99), 

au (0:99,0 : 99), p(0:99, 0:99), pf (0:99, 0:99), 

pet(0:99),pft(0:99),q(0:99,0:99), 

qext ( 0:99,0:99), 

qmax (0:99) ,rey(0:99,0:99), 

rho(0:99, 0:99), tamb, tfinal, time (0:99), 

u(0:99, 0:99), v(0:99, 0:99), 

vet ,x ( 0 :99 , 0 :99 ) , 2 ( 0:99) 


common 


6 


cp , ®t , f ram® , kc , nod® , pump , 
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a cpr,cvo, 

a d, deltat , do , • , aps , 

a h , he , hO , f , k , 

a k co , kel vin , ke 1 wa 1 , 

a l,m,mdot, 

a mu,p,pf,pet, 

i pft,q,qext,qmax, 

k ray , rho , tamb , tf inal , 

k tins , u , 

k v , vat , x , z 

if((frame .aq. 0) .or. 

& (node .aq. 0) .or. 

k (node .aq. at)) than 

if ( node .It. cp ) kal vin ( node , f rase ) = 20 . 39d0 
if(node .ge. cp ) kal vin ( node , f rame )atamb 
endif 

rho ( node , frame ) =p ( node , f ram a)/ 
k 41 57 . dO/kel vin ( node ( f raae ) 

u(node,frame)=2.5dO*4l57.dQ*kelvin(node,fraae) 
call propurho 
return 
end 

subroutine propurho 

integer*2 cp,et, frame, kc(5, 0:99), node, pump 
real * 8 cpr(0:99, 0:99), evo (0:99.0:99), 
k d(0:99),deltat, do (0:99), e{0:99, 0:99), eps<0:99), 

k h(0:99, 0:99), he (0:99, 0:99), 

k h0(0 :99,0 :99),f(0 : 99,0 : 99),k(0 : 99,0: 99), 

k kco (0:99, 0:99), kelvin (0:99, 0:99) , kelwal ( 0 : 99 , 

a 0:99), 1(0:99), m(0:99, 0:99) , adot ( 0 : 9 9 , 0 : 9 9 ) , 

k mu (0:99, 0:99), p(0:99, 0:99), pf (0:99, 0:99), 

4 pet ( 0 : 99 ) , pf t ( 0 : 99 ) ,q(0:99,0:99) , 

a qext ( 0: 99 , 0: 99 ) , 

a qmax ( 0 : 99 ) , rey ( 0 : 99 , 0 : 99 ) , 

a rho (0:99,0:99 ) , ta mb, tfinal, time (0:99) , 

a u(0:99, 0:99), v(0:99, 0:99), 

a vet , x ( 0 : 99 , 0 : 99 ) , z ( 0 : 99 ) 

common 

a cp , at , f rame , kc .node , pump , 

a cpr.cvo, 

a d, deltat , do , e , aps , 

a h ,hc , hO , f , k , 

a kco , kel vin , kelwal , 

a 1 , m , mdot , 

a mu , p , pf , pet , 

a pf t ,q ,qext , qmax , 

a rey , rho , tanb , tf inal , 

a time.u, 

a v, vet , x , z 

kelvin ( node., frame ) =u (node , f rame )/2 . 5d0/4157 . dO 
p ( node , frame )«rho (node, frame )*4157.d0* 
a kel vin ( node , f rame ) 

h(node, fra me )*1.4d0*u(node, frame) 
mu(node,frame)«2.d-6 
k(node,frama)=0.05d0 
return 
end 

block data h2block 

integer*2 cp,et, frame , kc ( 5 , 0 : 99 ) , node , pump 
real * 8 cpr(0:99,0:99), evo (0:99,0:99), 
a d(0:99), deltat, do (0:99), e(0:99,0:99) ,eps (0:99) , 

a h(0:99, 0:99), hc(0:99, 0:99), 
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a hO (0:99,0: 99), f(0:99, 0:99), k (0:99, 0:99), 

a kco( 0:99, 0:99 ) , kel vin < 0 : 99 , 0 : 99 ) , kalwal ( 0 : 99 , 

& 0:99), 1(0:99), n(0:99,0:99 ) , odot ( 0 : 99 , 0 : 99 ) , 

a au (0:99, 0:99), p<0: 99, 0:99), pf (0:99, 0:99), 

a pet (0:99 ) ,pft(0 :99) ,q(0:99,0 :99 ) , 

a qext (0:99,0:99), 

a qmax (0:99) , ray (0:99, 0:99) , 

k rho (0:99,0:99), tamb , tf inal , tioa ( 0 :99 ) , 

k u(0 :99 ,0:99 ) , v(0:99 ,0 :99 ) , 

a vet, x(0:99, 0:99) ,z(0:99) 

COBDOIt 

k cp , at .frame ,kc , node , puap, 

k cpr,cvo, 

k d, del tat , do , a , eps , 

k h ,hc , hO , f ,k , 

& kco,kelvin, kalwal, 

k 1 , 0 , odot, 

k au , p , pf , pet , 

k pft,q, qext, qmax, 

k ray , rho , tamb, tf inal , 

k time,u, 

a v, vat , x , z 

data cp , at/4, 35/ 
data d/ 

k 4*0. 2d0, 32*0. 25d0, 64*0. d0/ 

data do/ 

k 4*0. 22d0, 32*0. 27d0, 64*0. d0/ 

data eps/ 

& 100*0. 0d0/ 

data 1/ 

k 100 * 10 . do/ 

data qmax/ 

k '100*1500. d0/ 

data z/ 

k 10. dO, 34*0. dO, 30. dO, 64*0. d0/ 

end 


subroutine trans 

integer *2 cp , et , frame , kc ( 5 , 0 : 99 ) .node , puap 
real *8 cpr(0:99,0:99) , cvo ( 0 : 99 , 0 : 99 ) , 
k d(0:99),deltat, do (0:99), a(0:99, 0:99) , eps (0:99) , 

k h(0 :99, 0:99 ) ,hc(0:99, 0:99) , 

k hO (0 :99, 0:99), f(0:99, 0:99), k(0:99, 0:99), 

k kco( 0:99, 0:99) , kal vin ( 0 : 99 , 0 : 99 ) .kalwal ( 0 : 99 , 

k 0:99), 1(0:99), b( 0:99, 0:99) , adot ( 0 : 99 , 0 : 99 ) , 

a mu (0:99,0: 99), p(0:99, 0:99), pf (0:99, 0:99), 

a pet ( 0 : 99 ) , pf t ( 0 : 99 ) ,q(0:99,0:99) , 

a qext ( 0 : 9,9 , 0 : 99 ) , 

a qmax ( 0 : 9'9 ) , rey ( 0 :99 ,0 :99 ) , 

a rho ("0 :99, 0:99), tamb, tf inal , time (0:99), 

a U(0:99, 0 :99 ) ,v(0:99, 0:99) , 

a vet, x(0:99, 0:99) ,z(0:99) 

common 

a cp , et , f rame , kc , node ,pump , 

a cpr.cvo, 

a d , del tat , do , e , eps , 

a h , he ,h0 , f ,k , 

a kco , kal vin , kalwal , 

a l,m,Bdot, 

4 mu , p , pf , pet , 

a pft,q, qext, qmax. 
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k 

k 

k 


k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 


r ey , r ho , t amb , t f ina 1 , 
time,u, 
v , vat , x , z 

f(node,frame)=0.002d0 
he ( nods , 0 ) =0 . dO 
hc(node,frame)«1360.d0 

if ( node .It. cp ) kelwal ( node , 0 ) =20 . 39d0 
if(node .ge. cp ) kel wa 1 ( noda , 0 ) = t aab 
qext ( node , f rame ) =qmax ( node ) 

• ( tamb-kelwal (node , f rane-1 ) ) 

/ ( t a mb- 20 . 3 9d0 ) 

kelwal (node , frame )«>kelwal (noda , frame-1 ) 

+qext ( node , f rame-1 ) 

/4 19/7 8 00/da tan ( 1 . dO ) 

/( do ( node ) * * 2 
-d(node)**2)/l(node) 

* de 1 1 a t-hc ( node , f r aae-1 ) 

*4 . dO * da tan ( 1 . dO ) *d ( node ) *1 ( node ) 
* ( kelwal (node , frame-1 ) 
-kelvin(node, frame-1) ) 

•deltat 

/4 19. dO/78 00 .dO/da tan (l.dO) 

/ ( do ( node ) * *2 
-d(node)*‘2)/l(node) 


return 

end 

subroutine prophp 

integer‘2 cp, at, frame, kc(S, 0:99), node, pump 
real *8 cpr(0:99,0:99) , evo ( 0 : 99 , 0 : 99 ) , 
k d(0:99), deltat, do (0:99), 6(0:99,0:99) , eps (0:99) 

k h(0 :99, 0 :99 ) ,hc(0 :99, 0:99 ), 

i hO (0:99, 0:99), f(0:99, 0:99) ,k(0:99, 0:99), 

k kco( 0:99,0:99) , kel vin ( 0 : 99 , 0 : 99 ) , kelwal ( 0 : 99 , 

k 0:99), 1(0:99), m(0:99, 0:99) , adot ( 0 : 99 , 0 : 99 ) , 

k mu (0:99, 0:99), p(0:99, 0:99), pf(0: 99, 0:99), 

k pet(0:99 ) ,pft(0:99) ,q(0:99,0:99) , 

k qext (0:99,0:99), 

a qmax ( 0 : 99 ) , ray ( 0 : 99 , 0 : 99 ) , 

k rho ( 0 : 99 , 0:99), ta mb, t final , time (0:99), 

k u(0:99, 0:99), v(0: 99, 0:99), 

4 vet, x(0:99, 0:99) ,z(0:99) 


common 


k cp , at , f rame , kc , noda , pump , 

k cpr,cvo, 

k d , del ta t , do , e , eps , 

k h , he , hO , f , k , 

k kco , kel vin , kelwal , 

k l,m,mdot, 

k mu , p , pf , pet , 

k pf t , q , qext , qmax , 

k rey , rho , tamb , t final , 

k time , u , 

k v , vet , x , z 

u ( noda , frame )«h(node, frame )/l. 4 dO 
rho ( node , frame ) »p ( node , f rame ) / ( h ( node , frame ) 
k -u ( node , f rame ) ) 

call propurho 
return 
and 

subroutine feedtank 

integer*2 cp, at, frame, kc(5, 0:99), node, pump 
real*8 cpr(0:99,0:99),cvo (0:99, 0:99) , 
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k d<0:99), deltat, do (0:99), *(0:99,0:99) , eps (0:99) 

t h(0:99, 0:99), he (0:99, 0:99), 

k hO (0:99, 0:99), £(0:99, 0:99) ,k(0: 99 ,0:99), 

k kco( 0: 99,0: 99 ) ,kelvin( 0:99,0:99) .kelwal <0:99, 

k 0:99), 1(0:99), >(0:99,0:99) ,mdot ( 0 : 99 , 0 : 99 ) , 

k >u (0:99, 0:99), p(0: 99, 0:99), pf (0:99, 0:99), 

k pet (0:99) ,pft(0:99),q(0:99,0:99), 

k qext ( 0:99,0:99) , 

& qmax ( 0 : 99 ) , rey ( 0 : 99 , 0 :99 ) , 

k i ho (0 :99, 0 :99 ) , tamb , tf inal , tine ( 0 : 99 ) , 

k u(0 : 99, 0:99), v(0 -.99,0:99), 

k vet, x(0:99, 0:99) ,2(0:99) 

COBBOn 

k cp , et , f rame , kc , node , pump , 

k cpr, evo, 

k d , deltat , do , e , eps , 

k h , he , ho , f , k , 

k kco , kel vin , kelwal , 

k 1 ,m,mdot , 

k BU , p , pf , pet , 

k pf t ,q,qext , qmax , 

i rey , rho , tamb , t final , 

k time , u , 

k v , vet , x , z 

2 ( 0 )=z ( 0 ) 
return 
end 

subroutine extank 

integer*2 cp,et, frame, kc(5, 0:99) , node , pump 
teal *8 cpr (0:99,0:99) , evo ( 0 : 99 , 0 : 99 ) , 
k d(0:99), deltat, do (0:99), e(0:99, 0:99) , eps (0:99) 

& h(0:99, 0:99), he (0:99, 0:99), 

k hO (0:99, 0:99), f(0:99, 0:99), k(0:99, 0:99), 

k kco (0:99, 0:99) , kel vin ( 0 : 99 , 0 : 99 ) , kelwal ( 0 : 99 , 

«■ 0:99), 1(0:99), b( 0:99, 0:99) , mdot ( 0 : 99 , 0 : 99 ) , 

k mu (0:99, 0:99), p(0:99,0: 99), pf (0:99, 0:99), 

k pet ( 0 : 99 ) , pf t ( 0 : 99 ) ,q(0:99,0:99) , 

k qext (0:99,0:99), 

k qmax ( 0 : 99 ) , rey ( 0 : 99 , 0 : 99 ) , 

k rho (0:99,0:99), tamb , t final , time (0:99), 

k u(0 : 99 ,0 : 99 ) , v(0 :99 ,0 :99) , 

k vet, x( 0 :99, 0 :99 ) ,z(0:99) 

common 

a cp , et , -frame , kc , node , pump, 

k cpr , evo , 

k d , deltat , do , e , eps , 

k h ,hc , hO , f , k , 

k kco , kel vin , kelwal , 

k 1 , a , adot , 

S mu , p , pf , pet , 

S pf t , q , qext , qmax , 

k rey , rho , tamb , tf inal , 

<> time,u, 

k V , vet , X , 2 

z ( et ) =z ( et ) 

return 

end 
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